home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 351-375 / disk_351 / pdc / libsrc.lzh / LibSrc / Math / sqrt.c < prev    next >
C/C++ Source or Header  |  1990-04-07  |  6KB  |  179 lines
  1. /************************************************************************
  2.  *                                                                      *
  3.  *                              N O T I C E                             *
  4.  *                                                                      *
  5.  *                      Copyright Abandoned, 1987, Fred Fish            *
  6.  *                                                                      *
  7.  *      This previously copyrighted work has been placed into the       *
  8.  *      public domain by the author (Fred Fish) and may be freely used  *
  9.  *      for any purpose, private or commercial.  I would appreciate     *
  10.  *      it, as a courtesy, if this notice is left in all copies and     *
  11.  *      derivative works.  Thank you, and enjoy...                      *
  12.  *                                                                      *
  13.  *      The author makes no warranty of any kind with respect to this   *
  14.  *      product and explicitly disclaims any implied warranties of      *
  15.  *      merchantability or fitness for any particular purpose.          *
  16.  *                                                                      *
  17.  ************************************************************************
  18.  */
  19.  
  20. /*
  21.  *  FUNCTION
  22.  *
  23.  *      sqrt   double precision square root
  24.  *
  25.  *  KEY WORDS
  26.  *
  27.  *      sqrt
  28.  *      machine independent routines
  29.  *      math libraries
  30.  *
  31.  *  DESCRIPTION
  32.  *
  33.  *      Returns double precision square root of double precision
  34.  *      floating point argument.
  35.  *
  36.  *  USAGE
  37.  *
  38.  *      double sqrt (x)
  39.  *      double x;
  40.  *
  41.  *  REFERENCES
  42.  *
  43.  *      Fortran IV-PLUS user's guide, Digital Equipment Corp. pp B-7.
  44.  *
  45.  *      Computer Approximations, J.F. Hart et al, John Wiley & Sons,
  46.  *      1968, pp. 89-96.
  47.  *
  48.  *  RESTRICTIONS
  49.  *
  50.  *      The relative error is 10**(-30.1) after three applications
  51.  *      of Heron's iteration for the square root.
  52.  *
  53.  *      However, this assumes exact arithmetic in the iterations
  54.  *      and initial approximation.  Additional errors may occur
  55.  *      due to truncation, rounding, or machine precision limits.
  56.  *      
  57.  *  PROGRAMMER
  58.  *
  59.  *      Fred Fish
  60.  *
  61.  *  INTERNALS
  62.  *
  63.  *      Computes square root by:
  64.  *
  65.  *        (1)   Range reduction of argument to [0.5,1.0]
  66.  *              by application of identity:
  67.  *
  68.  *              sqrt(x)  =  2**(k/2) * sqrt(x * 2**(-k))
  69.  *
  70.  *              k is the exponent when x is written as
  71.  *              a mantissa times a power of 2 (m * 2**k).
  72.  *              It is assumed that the mantissa is
  73.  *              already normalized (0.5 =< m < 1.0).
  74.  *
  75.  *        (2)   An approximation to sqrt(m) is obtained
  76.  *              from:
  77.  *
  78.  *              u = sqrt(m) = (P0 + P1*m) / (Q0 + Q1*m)
  79.  *
  80.  *              P0 = 0.594604482
  81.  *              P1 = 2.54164041
  82.  *              Q0 = 2.13725758
  83.  *              Q1 = 1.0
  84.  *
  85.  *              (coefficients from HART table #350 pg 193)
  86.  *
  87.  *        (3)   Three applications of Heron's iteration are
  88.  *              performed using:
  89.  *
  90.  *              y[n+1] = 0.5 * (y[n] + (m/y[n]))
  91.  *
  92.  *              where y[0] = u = approx sqrt(m)
  93.  *
  94.  *        (4)   If the value of k was odd then y is either
  95.  *              multiplied by the square root of two or
  96.  *              divided by the square root of two for k positive
  97.  *              or negative respectively.  This rescales y
  98.  *              by multiplying by 2**frac(k/2).
  99.  *
  100.  *        (5)   Finally, y is rescaled by int(k/2) which
  101.  *              is equivalent to multiplication by 2**int(k/2).
  102.  *
  103.  *              The result of steps 4 and 5 is that the value
  104.  *              of y between 0.5 and 1.0 has been rescaled by
  105.  *              2**(k/2) which removes the original rescaling
  106.  *              done prior to finding the mantissa square root.
  107.  *
  108.  *  NOTES
  109.  *
  110.  *      The Convergent Technologies compiler optimizes division
  111.  *      by powers of two to "arithmetic shift right" instructions.
  112.  *      This is ok for positive integers but gives different
  113.  *      results than other C compilers for negative integers.
  114.  *      For example, (-1)/2 is -1 on the CT box, 0 on every other
  115.  *      machine I tried.
  116.  *
  117.  */
  118.  
  119. #include <stdio.h>
  120. #include "pml.h"
  121.  
  122. #define P0 0.594604482                  /* Approximation coeff  */
  123. #define P1 2.54164041                   /* Approximation coeff  */
  124. #define Q0 2.13725758                   /* Approximation coeff  */
  125. #define Q1 1.0                          /* Approximation coeff  */
  126.  
  127. #define ITERATIONS 3                    /* Number of iterations */
  128.  
  129. static char funcname[] = "sqrt";
  130.  
  131. extern double frexp ();
  132. extern double ldexp ();
  133.  
  134. double 
  135. sqrt (x)
  136.         double x;
  137. {
  138.     int k;
  139.     int bugfix;
  140.     int kmod2;
  141.     int count;
  142.     int exponent;
  143.     double m;
  144.     double u;
  145.     double y;
  146.     double rtnval;
  147.     struct exception xcpt;
  148.  
  149.     DBUG_ENTER ("sqrt");
  150.     DBUG_3 ("sqrtin", "arg %le", x);
  151.     if (x == 0.0) {
  152.         rtnval = 0.0;
  153.     } else if (x < 0.0) {
  154.         xcpt.type = DOMAIN;
  155.         xcpt.name = funcname;
  156.         xcpt.arg1 = x;
  157.         if (!matherr (&xcpt)) {
  158.             fprintf (stderr, "%s: DOMAIN error\n", funcname);
  159.             errno = EDOM;
  160.             xcpt.retval = 0.0;
  161.         }
  162.     } else {
  163.         m = frexp (x, &k);
  164.         u = (P0 + (P1 * m)) / (Q0 + (Q1 * m));
  165.         for (count = 0, y = u; count < ITERATIONS; count++) {
  166.             y = 0.5 * (y + (m / y));
  167.         }
  168.         if ((kmod2 = (k % 2)) < 0) {
  169.             y /= SQRT2;
  170.         } else if (kmod2 > 0) {
  171.             y *= SQRT2;
  172.         }
  173.         bugfix = 2;
  174.         xcpt.retval = ldexp (y, k/bugfix);
  175.     }
  176.     DBUG_3 ("sqrtout", "result %le", xcpt.retval);
  177.     DBUG_RETURN (xcpt.retval);
  178. }
  179.